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АВ5ТКАСТ 

The general solution of the behavior of a viscously damped vibrating 
system having several degrees of freedom, as developed in the literature 
via matrix methods, is treated in detail here and made the basis of a 
digital computer program which is capable of determining the natural 
frequencies, the mode shapes, and the displacements of each mass as 
functions of time. This program, which is written in FORTRAN 60 language 
for the Control Data Corporation 1604 computer, affords several output 
options. It does not treat cases of supercritical damping or cases in 


which two or more natural frequencies are the same. 
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CHAPTER 1 

INTRODUCTION 
1.1 General Remarks. The analysis of a subcritically damped multiple 
degree of freedom vibrating system necessitates obtaining the solution of 
a complex eigenvalue problem to determine the natural frequencies and mode 
shapes of the system. Although the analysis presented in the literature 
for systems with two degrees of freedom may be extended to systems with 
more than two degrees of freedom, manual calculations are too laborious 
to be practical. Therefore the natural frequencies are usually found by 
ignoring the presence of the damping and solving the resulting real eigen- 
value problem. This simplifies the problem considerably and provides a 
good approximation provided the damping is light. Another technique 
employed is to solve the real eigenvalue problem and then account for 
the damping using perturbation techniques. However, even in the absence 
of damping, systems involving more than two degrees of freedom usually 
require iteration or trial and error techniques (such as the Stodola or 
Holzer methods) to obtain the mode shapes. 

The advent of the electronic digital computer has eliminated the 
necessity of ignoring the damping component and increased the size of the 
system for which solutions can be obtained. Although not entirely devoid 
of error the digital computer is highly reliable and its speed of opera- 


tion has made it an invaluable tool in engineering analysis and design. 





1.2 Method of Analysis. The multiple degree of freedom damped vibrating 
system is described by a set of N linear second order differential equa- 
tions, where N denotes the degrees of freedom involved. Utilizing matrix 
analysis and generalized coordinates the system is described by 

МХ + СХ + KX = FQ) 
which is the same form as the single degree of freedom system. However, 
M, C, and K now represent square matrices and X, X, X, and F(t) are 
column vectors. The analysis presented in what follows is performed using 
matrix technique because of the compactness of the notation and the order- 
ed computational procedure which is ideally suited for digital computer 
programming. 

The mathematical analysis is demonstrated and equations are derived 
which describe the time behavior of free and forced vibrating systems. 
Finally a digital computer program is presented which performs the opera- 
tions indicated in the mathematical analysis to determine the natural 


frequencies, mode shapes, and time behavior of the damped vibrating system. 





СНАРТЕК 11 
MATHEMATICAL ANALYSIS 
2.1 Eigenvalues and Eigenvectors. The viscously damped vibrating system 
with one degree of freedom is described by a linear, second order differ- 
ential equation. 
MK +CX + KX= F(t) 
M, C, and K represent the mass, damping and stiffness of the system, and 
F(t) the forcing function. A system with many degrees of freedom is de- 
scribed by a set of second order differential equations similar to the 
case of the single degree of freedom system. Using matrix notation the 


system is described for the case of free vibration by 
(1) MX + CX + KX =O 


where M, C, K, are now square matrices of order N, representing, as in 
the single degree case, the mass, damping and stiffness matrices of the 
system. It is always possible to write the equations so that these mat- 
rices are symmetric. М 1з the number of degrees of freedom. X, X, and X 
are column vectors of order N representing displacement, velocity, and 
acceleration in generalized coordinates. 

Premultiply equation (1) by (en letters will hereafter re- 


present matrices and lower case letters scalar constants) to obtain; 
(1а) X + MCX+ MKX =O 


* 
W. J. Duncan and A. R. Collar [1] have shown a simplified method of find- 
ing the eigenvalues and eigenvectors by writing the second order differen- 
tial equation in reduced form as a first order differential equation. The 
m 


Numbers in brackets refer to references listed in the bibliography 
page 23 





reduced form is obtained by introducing velocity as a dependent variable. 


Let 
X 
(2a) с B 


Here Y is a partioned column matrix of order 2N with the N displacement 
components in the upper half and N velocity components in the lower half. 


Let 


(2b) je | > c 
-M'K -м'с 


Here 0 is an nth-order null matrix and I is an nth-order identity matrix. 


The reduced equation then becomes 
(3) Y= AY 


A solution of a first order differential equation such as equation 


(3) may be taken as 


(4) MES CE 


Substitution of the assumed solution into equation (3) gives 


τ et 


(5) pe U = Ae u 


Dividing out the exponential term and rearranging, we have 
(5a) [а = O 


To have a nontrivial solution, it then follows that 

det lA - pT] =O 
This is the characteristic matrix of the reduced system. The eigen- 
values or roots of the characteristic equation are Pio = 1; 2; 5 2N 


and the eigenvectors are U; т = νο... Ν 





where 
Uir 


) 
Mar 
Up = \ 


Uan] 


The first subscript indicates the position in the column and the second 
subscript corresponds to the rth eigenvalue. 

As long as the damping is less than critical in each mode, the re- 
duced system will yield N complex conjugate pairs of eigenvalues and 
eigenvectors. The eigenvalues will be of the same form as the sub- 
critically damped single degree of freedom system. 

Pr = - u. + j Or 1-52 
Pr = — „оо, - jryi- SF 
The bar indicates the complex conjugate and the subscript indicates the 


mode, and where 


damping ratio 


5, 


ملف 


natural frequency 
a, Y1- - damped natural frequency 
5.08 - decrement factor 

The corresponding eigenvectors are 

ur = M +0 Wr 

Ur = Vr = j We 
2.2 Homogeneous Equation. 
221 Orthogonality Relations. In the case of the undamped vibrating 


system it is always possible to choose a set of coordinates in which the 
mass and stiffness matrices are uncoupled. However in the damped vibrat- 
ing system, unless the damping matrix is proportional to either the mass 

or stiffness matrix, a set of coordinates which will uncouple the equations 


5 





of motion cannot be found without a knowledge of the eigenvalues and eigen- 
vectors of the system. 

K. A. Foss has developed a set of orthogonality relations for the 
eigenvectors of a damped linear dynamic system from which coordinates 
which uncouple the equations of motion may be found. [2] The technique of 
reducing the second order differential equation to one of the first order 
is again employed, with slightly different notation. 


Let 


о м ІС : 
(6a) .-|9 " (6b) TP 4 (Come K 


Using the above notation, equation (1) becomes 


7) RF + BZ =O 


Furthermore let 


Pr 
X = Ure 
| e. t 
X = р. Ч; ё 
Then 
p.t 
(8) z = me 
where = pror 
OF Sie 


Substituting equation (8) into equation (7) and dividing out the exponen- 


tial factor we get, 


(9) pr RO, + B§, =0 


Since the eigenvalues and eigenvectors occur in complex conjugate pairs 


the complex conjugates of Us and ps also satisfy the homogeneous equation. 





Denoting the complex conjugate by the subscript S, equation (9) may be 


written 


(10) psRO, + Bos =o 


Or since M, C, and K are symmetric, R and B must be symmetric, and equa- 


tion (10) may be written 
т T 
(10a) Ps de R + $. B = О 


T 
Premultiply equation (9) by Ф. , postmultiply equation (10a) by Ф, 


and subtract (10a) from (9). 
T 
p, &R $, « 8,89,=0 
- e, & RS, + ВФ, = о 


(Р - рз) PRÈ, gka 





Therefore unless r = s the orthogonality relations are 
т 
11) $,RÓ$, -ο 


an 818%, = о 


20212. General solution of the homogeneous equation. In the general 
solution of equation (7), 2N arbitrary constants of integration must be 
evaluated. The 2N initial conditions of displacement and velocity are 
used to evaluate these constants. Assume a solution of the form 
ZN 
pt 
a3) z(= > $ C, e 


ға! 
where C. represents the 2N constants to be determined. 
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The initial conditions may be expressed as a vector expansion of the 


natural modes 


2N 
(14) z(-) qu 
ТТ 


т 
Premultiply both sides of equation (14) by Ф. R 


ZN 
(15) MITO Е >. $ RO, Cr 


Using the orthogonality relation, equation (11), the summation simplifies 


and equation (15) becomes 
1  Q RZ(Q9s PRÓ$,C 
(16) - = ртт 


Or 


un Ce = ғғ 
$ R $, 


The general solution is therefore found by substituting the value of the 


constants of integration into equation (13). 
ZN 


2 


It follows from Ee evaluation of equation (18) at time zero, that 


ZN T 
Y GER. 1 
= $- R $, 
Equation (18) may be written in the following partitioned form using the 


nth-order system notation. 





, = ὑπο 2 r P, U-ULMX(d + πε X (0) 
(19) [ 


B о 
XQ) e UrULMX6) + unur Cx6) 


Ml 
ὁ am ie + UL CU, 


ZN 


e.t 
Е 


from which 


NC luu Mx Co) + Pr Ur Ur MX Co) + П) ert 


(20) a: 2р. ці ми. + ul Cup 


Dei 
8 





It should be noted again at this point that for the sub-critícally damped 
system being considered, pu and U are complex quantities and occur in 
complex conjugate pairs. The matrix operations indicated in equation (20) 
have been carried out by breaking up the complex quantities into real and 
imaginary parts and then recombining. The operations are relatively 


straightforward but lengthy, therefore only the final form will be pre- 


sented. 
N 


2 3 
)21( x» | eta + Cam, Н.-М +6.С)хсо) Cos в.) е 


Pel 


N 
): ΄ 8 
2 Я = Ap 
+ Eu μενῶ (е. Нм -B,GrM E Wo X e) — © 


еті 


Where: 
Prz ¬ Kr t jêr U. * Ve t j We 


αγξ - Κας (WMV. 3 Wr MWr) = #Pr V MW. + у; СУ, т We C Wr 


+ E 
b= 2 Be CVEMVe - WI Mw) - 4% Ve Mi + £ Ve C Wr 


Ae VENE - WWE 
By = VW +WrYre 
Gr = GrAr +byB, 
H: = b- A, ó arBr 


The factor of two in the above equation results from the reduction of the 
summation from 2N to N. When the indicated operations of equation (20) 
are carried out from r = 1 to r = 2N the imaginary components sum to zero; 
therefore all quantities in equation (21) are real. Equation (21) has 
been given by S. H. Crandall and R. B. McCalley in reference 3. 
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2.3 Non-Homogeneous Equation 
2-3-T Steady State Solution of the Non-Homogeneous Equation. Once 


again using the reduced form, the non-homogeneous equation becomes 


(22) RZ +Bz = F(>) 


© 
Fs FA 


In a linear nth-order differential equation the solution of the non- 


Where 


homogeneous equation may be found, once the solution of the homogeneous 


equation is known, by expanding 2 in a modal series. (2,4] 


Thus, ZN 
(23) Е(д- D %,У- 0) 
rs | 
here  Y.() is a variable parameter. 


Upon ασε εσείς of equation (23) in equation (22), we obtain 
2. 


(24) > REY, + Bê = FOA 
ES T 
Premultiplication of equation (24) by p and making use of the 


orthogonality relations, yields 
m ë T т 
(25) OROY+ $ BEY = $ FO 


Or 


(25a) Rn Y = Pr Kn Y = БА 
Where 
E 
Ra = $ R $, 
=PrRn = $. 86. 
m= $ FQ) 
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Using the Laplace Transform technique,the solution of equation (25a) 


is easily determined. 





Rearranging 
F. 
(a) Y - p. Y = T 
The Laplace Transform is 
Fa CS 
(b) sY($) - pr Y CS) Ry 


Therefore 


| ҒА (6) 
(с) ҮСӘ) = R a 


The inverse transform of (c) is 


+ 
- (¥-%) 
(26) YC) EL F. @ d< 


Since pr = - Xr + Pr 
t 
- potet 
(264) Y(t) = ЗЕ š da (^c) (cos e. (+-©) +) 5м 6 (1-)) 4 τ 


Substitute (26a) into (23) 


Z t 
W a 
(27) zQ) NES | Έω (οο»β, (1-1) +5 sin (сч 
r=! 2 


As in the case of the homogeneous solution the reduced system quantities 
are replaced by their nth-order equivalence. The final form of the 


steady state solution then becomes 


t 
-Xe t-a) 
(28) ТЭР [εως a Cos 4, (t-T) AT 


2 f f) e r SIN By (%-2 AT 


The notation in equation (28) is the same as that used in the homogeneous 





solution. 
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2.902 Steady State Solution with Sinusoidal Excitation. A special, 


but very common problem is the case in which the forcing function, { (9 А 


is sinusoidal.  Expressing the forcing function as 
f (сай) 
(29) (y so (De 


where @ - indicates "the real part of" 


D - column of driving force amplitudes, possibly complex to admit 
phasing 


0% - excitation frequency 


The non-homogeneous equation therefore becomes 
e А ГАЙ 

(30) МУ + СХ +Kx = Ф (ое 

Assume a solution of the form 


ТЕ. 
G1) x) = @ (Ge) 


where G із а column of undetermined constant coefficients. Substitute 
the assumed solution in equation (30) and divide through by the exponen- 


tial factor to obtain 
(32) -W MG + j CG +KG =D 
Or 
=] 
(33) G= [к-ам 25 


Therefore the steady state solution for the case of sinusoidal excitation 


is 


2 
ος 
(34) x (+) = elfk- am + CD el^ ] 
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2.4 Transient Plus Steady State Solution. 
2.4.1. Initial Conditions. The total solution is given by 


(35) XÐ = TS +SS 


where TS is the so called transient solution and SS the so called steady 
state solution. However in the total solution the initial conditions 

of displacement and velocity in the transient solution must account for 
the initial values in the steady state solution. This is accomplished 
as follows. In the case of sinusoidal excitation the steady state 


solution was found to be 
-l ‘ot 
(34) X (+) = е! [к-ам +juwC|D e | 


To allow for phase differences in the exciting forces on the masses, let 


D be given by 

(36) D = DR + ¡DI 

and let the inverse of |κ- wM +j wC] be denoted by 
/ | / 

(37) A +j 


jet 
Then since € = Coswt + ; sınwt, equation (34) 


becomes 
(38) xt = 2 [+53] (D+; o1) ( ces ust +) SIN a?) 


Carrying out the indicated multiplications, we obtain 


(39) X) = [ (on - oz) cos «t - (BDR + ADI) SIN ωχ] 
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and the velocity is given by 
(40) X (*%) = [- t» (ADR-B'DI) sINWt — wW(BDR + A’DT) Cos w қ 


Therefore the initial values of displacement and velocity of the steady 


state solution are obtained by evaluating equations (38) and (39) at t=o 


(41) X. = ADR-BDI 


(42) X0) = -w(BDR+ADT) 


The initial conditions used in determining the time behavior of the 


transient portion for the total solution therefore become 


(43) xX, @) = Xx) -(A'DR-—BD1) 


" 


(&&) Y, (о) 


і! 


X (o) + о (втв + АРТ) 


The complete time-solution is thus obtained as indicated in equation (35) 
where the transient solution is found by using the modified initial condi- 


tions of equations (43) and (44). 
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CHAPTER III 
PROGRAM DESCRIPTION 

3.1 General Remarks. A digital computer program "PROGRAM VIBSYS" is 
presented which performs the mathematical operations of the equations 
developed in the previous chapter. The program is coded in Fortran 60 
programming language, [s, 6] specifically for the Control Data Corpora- 
tion 1604 computer. Although the Fortran 60 language is applicable to 
most large digital computers there are minor variations which are pecul- 
iar to the specific system in use. The Fortran 60 language does not per- 
mit automatic operations with complex numbers; therefore all operations 
involving complex numbers are accomplished by operating on the real and 
imaginary parts separately and then recombining. 

The eigenvalues and eigenvectors of the reduced system are found by 
a matrix iteration scheme which utilizes the direct and inverse power 
methods and matrix deflation. A mathematical subroutine MATSUB is used 
to carry out these operations. (7,8) Ав presented, a maximum of twenty 
iterations will be performed using the direct power method and then a 
maximum of twenty using the inverse power method. If convergence is not 
reached with the inverse power method a print out to this effect will be 
executed. 

Subroutine "INVERT" is used for matrix ene tame "INVERT" uses 
the Gaussian elimination and pivotal techniques. Inversion of the 
[<- a M+ ©] matrix which contains complex elements is achieved in 

* 


Subroutine INVERT is a library subroutine of the Naval Postgraduate 
School and is designated locally as Fl NPS INVERT, 
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the following manner. Let n + jB) be the matrix to be inverted and 
[ο + 3D] the inverse to be determined. Then by definition 
La +58] [c +0] = І 
The above matrix equation leads to two simultaneous equations with 
two unknowns, C and D. 
(a) АС-Вр = І 
(b) BC+AD= © 
Solving first for C 


C 


А Кп 
ΠΠ 


and then for D 


- ABC 


D 

Thus the complex inversion is reduced to multiplications, additions, 
and inversions, of matrices of real numbers. 

Flexibility and utility are the principal aims of the program. 
Usage requires a knowledge of the mass, stiffness, and damping matrices. 
Although various output options are available which require additional 
input data, the three above mentioned matrices are all that are necessary 
to determine the eigenvalues and eigenvectors of the reduced system and 


the natural frequencies and mode shapes of the original system. 


3.2 Program Options. In addition to finding the natural frequencies and 
mode shapes of the system, five options are available which describe the 
time behavior of the system under conditions of free and forced vibration. 
(a) Option l. The time solution of the free vibration problem 
is obtained in general form and no additional input data is required for 
execution of this option. In section 2.2.2 the general solution of the 


homogeneous equation was developed and is repeated here for convenience. 
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(21) XC) (в МХ (о) + (- 1 G-M+ BHM + Gr σα "gt 
a 
) La 2497 (кә БЕ <. Н.М - 4. &6-М + Ἡ 26) xe) |e E SIN prt 
rz 
The output of option 1 consists of the four coefficient matrices of the 
ха» со5 Ж; X Co) COS Bet ) X (o) SIN brt and `X Co) Sin Be t 
terms. Therefore, the output of option 1 will consist of 4N square matrices. 
(b) Option 2. The execution of option 2 requires as additional 
input data the values of the initial displacement, and initial velocity 
vectors. The product of the coefficient matrices of option 1 and the 
initial displacement vector or initial velocity vector, as appropriate, is 
performed to obtain a coefficient column for the cosine and sine terms. 
(c) Option 4.* The general steady state solution of the forced 
vibration problem is provided by option 4. The general solution with a 
forcing function f(t) was developed in section 2.31 and found to be 


N 


* 
= S (t-t 
(28) x(t)= 27 feme Ë cos С®-*®) ат 


"sl " 


5 
o) die f fe SIN Gr (r-*)act 
ос + ©, 2 


ral 
In this case the coefficient matrices of the convolution integral are 


2 


evaluated. There will be 2N square matrices, N corresponding to the 
coefficient of the integral involving the cosine term and N coefficient 
matrices of the integral involving the sine term. 
(d) Option 5. The steady state solution of the special case 
of forced vibration with sinusoidal excitation is determined in option 5. 
* 
Option 3 is described in subsection (e), following. 
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In section 2.4.1 the expression for the time behavior was found to be 


4 


(39) X() = (KoR- Bor) coswt - (BDR+ADT) sınwt 


Additional input data required for the execution of this option include 
the driving force amplitude and excitation frequency. The output consists 
of the coefficient column vectors of the cosine and sine terms. 
(e) Option 3. A plot of displacement versus time is made for 
each mass, for one of the three following cases. 
1. Free Vibration 
2, Steady state vibration with sinusoidal excitation 
3. Transient plus steady state with sinusoidal excita- 
tion 
Case l is obtained when the initial conditions of displacement and veloc- 
ity are given as input data and the driving force amplitudes and excita- 
tion frequency are zero. 
Case 2 is obtained when the driving force amplitudes and excitation 
frequency are given and the initial conditions of displacement and 
velocity are zero. 
Case 3 is obtained when the initial conditions, driving force amplitudes 
and excitation frequency are all given. 
Any combination of the options may be obtained with one set of input 
data. Input data format and output control is described in detail in 


Appendix B. 


3.3 Accuracy of Method. The accuracy of the solution is contingent on 
the accuracy with which the eigenvalues and eigenvectors are determined 
and the effect of computer roundoff error. Internal checks have been 


provided so that the integrity of the solution may readily be evaluated. 
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3.3.1 Internal Checks. Prior to the determination of the eigenvalues 
and eigenvectors the trace and determinate of the matrix, A, is calcu- 
lated. A cia therefore readily available since the trace of the 
matrix A is equal to the sum of the eigenvalues of the characteristic 
matrix of A and the determinant is equal to the product of the eigen- 
values. The trace and determinant of A and the sum and product of the 
eigenvalues are included as standard output in Program VIBSYS. 

Additional checks are inherent in the solution. In Option 1 the 
summation of the coefficient matrices of the Xx(9 Cos,S-t terms 
must sum to the identity matrix, since the evaluation of x(t) at time 
zero must equal the initial conditions of displacement. This is easily 
seen by referring to equation 21 and considering the initial values of 
velocity to be zero. Similarly, in Option 2 the sum of the coefficient 
columns of the cosine terms must equal the initial conditions of dis- 
placement. 

In the steady state solution with sinusoidal excitation,the inversion 
of the complex matrix presents one of the greatest possible sources of 
error. Therefore when Option 5 is executed, the output includes the real 
and imaginary parts of the inverse and the product of the original matrix 


and its inverse. This product should, of course, be the identity matrix. 


3.3.2 External Checks. In order to determine the reliability of 
Program VIBSYS, sample problems of two degrees of freedom for which solu- 
tions were available in the literature were programmed. Correlation of 
results was excellent. No suitable sample problems involving more than 
two degrees of freedom were found in the literature. Therefore, in order 
to extend the range of testing, solutions of undamped systems were compar- 


ed with program solutions of similar systems with negligible damping. 
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The results were as anticipated; as the damping was decreased the natural 
frequencies approached those of the undamped system. The results of a 


sample problem are demonstrated in Appendix E. 


3.4 Program Limitations. Although the digital computer extends the size 
of system for which solutions are obtainable, the ultimate size is dictate 
ed by the storage capacity of the computer. The CDC 1604 computer has 
approximately 32,000 storage locations; however, it is easily seen that 
this may be rapidly exhausted. Program VIBSYS is limited to a system 
with 10 degrees of freedom. 

The basic principle in the development of the general equations is 
the orthogonality relations of the eigenvectors and the parameters of 
the system. Therefore each eigenvalue must have associated with it a 
unique eigenvector. This requirement is not fulfilled in the case of 
a system having two equal eigenvalues, and therefore two natural fre- 
quencies of equal magnitude. In such a case VIBSYS determines the 
frequencies and the program terminates at that point. 

If the damping is critical in any part of the system an attempt 
will be made to obtain the eigenvalues and eigenvectors. However, if 
they are found, they will not occur in the complex conjugate pairs and 


the program will be terminated. 
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CHAPTER IV 
CONCLUSIONS AND RECOMMENDATIONS 

Program VIBSYS provides an accurate method for determining the natural 
frequencies, mode shapes, and time behavior of a subcritically damped 
dynamic system. The compiling time for the program is 4 minutes and 30 
seconds, while the running time for a system with four degrees of freedom 
is approximately one mínute. (Running time will depend on the number of 
iterations required to obtain the eigenvalues and eigenvectors.)  There- 
fore it is more efficient, with respect to computer time, to make multiple 
runs. 

Comparison of natural frequencies of undamped and lightly damped 
systems has shown that the undamped approximation is very good for damping 
ratios of less than 0.01. For systems tested, up to four degrees of 
freedom, the natural frequencies of the damped and undamped systems did 
not differ in the first three significant figures. 

The program may be extended in a variety of possible ways. The 
general nature of Program VIBSYS requires liberal use of computer stor- 
age space. The size of the system may be extended by segmenting the 
program into separate programs for handling specific problems such as, 
free vibration, forced vibration with a general forcing function, or 
forced vibration with sinusoidal excitation. Furthermore, with minor 
modifications the present program may be terminated upon finding the 
natural frequencies and mode shapes and used as a separate program to 
obtain input data for specialized programs. 

Program VIBSYS can be augmented and made more useful by eliminating 


the restrictions which cause the program to be terminated if either (a) 
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the system has two eigenvalues of equal magnitude, or (b) one of the modes 


has aperiodic free motion, that is a purely real eigenvalue. 
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APPENDIX A 
PROGRAM STRUCTURE 

A.l General Remarks. Program VIBSYS is composed of a main body which is 
divided into five sections, and two subroutines, INVERT and MATSUB.  Func- 
tion subroutines used include the sine, cosine, exponential, square root, 
and absolute value functions. These function subroutines are called from 
the Fortran library tape. Subroutine draw is also available on the For- 
tran library type and is programmed specifically for the CDC 1604 computer 
for use with the CALCOMP plotter. 

Subroutine MATSUB, which is available as a CO-OP Library mathematical 
subroutine, required minor modifications in order to retain the eigen- 
values and eigenvectors in storage for use in the main body of the program. 
The flow chart provided with MATSUB in the CO-OP Library was found to be 
inadequate for understanding. In order to perform the necessary modifica- 


tions a revised flow chart was made up and included in Appendix C. 


A.2 Main Body. The main body consists of five sections with functions 
as listed below 

(a) Input - In addition to allocation of storage spaces for 
dimensioned arrays and setting up a control for reading input data, the 
input section contains a control sequence for MATSUB and curve labeling. 

(b) Natural frequencies and mode shapes - The reduced character- 
istic matrix is constructed and the eigenvalues and eigenvectors are deter- 
mined. The natural frequencies and mode shapes are then found. 

(c) Forced Vibration - The coefficient column vectors of the 
sine and cosine terms of option 5 are calculated. 

(d) Free Vibration - The outputs of options 1, 2, and 3 are 


formulated. 








-—— q^ Y) SD coma 
σαν... c à 
— СД 


ab мне 9 ων. 


(e) Forced Vibration with General Excitation - The coefficient 


matrices of option 4 are calculated. 


A.3 Subroutines. 
(a) MATSUB - Evaluates the eigenvalues and eigenvectors of the 
reduced system. 


(b) INVERT - Inverts a real matrix. 





ALIS 


ALRS 


BRN 


CD 


CDC 


CPD 
DI 
DISPI 
DISPR 
DM 
DR 
EP1 


EP2 


PROGRAM NOTATION 
coefficient column vector of cos P, t term, output of 
option 2 
coefficient column vector of cos4p t term when Х(0), апа 
X(0) are modified to obtain total solution 
control parameter for MATSUB 
control parameter for MATSUB 
inverse of mass matrix 
product of inverse of mass matrix and damping matrix 
mass matrix lb-sec?/in 
product of inverse of mass matrix and stiffness matrix 
-2 X (V, MV, = W TMW) -4 @, | v Tw, + v Tev, = w Tew 
imaginary part of identity matrix 
2 Êr (V TMV, - W Ta.) - 4 o v, Tu + 2v Tov, 
coefficient column vector of sin „© term, output of 
option 2 
coefficient column vector of sin @,t term when X(0) and 
X(0) are modified to obtain total solution 
coefficient column of sinwt, output of option 5 
imaginary part of driving force amplitude 
imaginary part of modal matrix A 
real part of modal matrix | 
damping matrix lb-sec/in 
real part of driving force amplitude 
iteration parameter for MATSUB 


iteration parameter for MATSUB 
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GBI 


GBR 


GRM 


HRM 


IC 


IDET 


IEG 


IFC 


ITITLE 


IVEC 


MIT 


MITS 


МР1 


МР2 


MP3 


MP4 


MP5 


RIN 


RPI 


SPEC 


STEP 


TIPI 


iteration parameter for MATSUB 

iteration parameter for MATSUB 

coefficient matrix of X(0) сов үс, output of option 1 
coefficient matrix of X(0) sin /S,.t, output of option 1 
control parameter for VIBSYS 

control parameter for MATSUB 

control parameter for MATSUB 

control parameter for VIBSYS 

title for graphical output 

control parameter for MATSUB 

controls number of iterations in MATSUB, power method 
controls number of iterations in MATSUB inverse power 
method 

option 1 control 

option 2 control 

option 3 control 

option 4 control 

option 5 control 

order of system 

real part of identity matrix 

coefficient column vector of cos c t term, output of 
option 5 

real part of the inverse of [к- a M PUJ w C) 
stiffness matrix lb/in 

spectral matrix 

step size for graphical output 

imaginary part of inverse of Ke ›2м + joc 
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UI 


VALI 


VALR 


VCO 


VCV 


VCW 


VECI 


VECR 


VO 


WCW 


WDM 


WO 


XCO 


XO 


хос 


XOS 


Y 


The above 


matrix of reduced system, imaginary part 

matrix of reduced system, real part 

imaginary part of eigenvalue of reduced system (4) 
real part of eigenvalue of reduced system (οι, 
initial velocity vector modified to account for initial 
velocity of steady state solution 
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Т 
Vp CW, 


V 


imaginary part of eigenvector of reduced system СУ) 

real part of eigenvector of reduced system (W.) 

initial velocity vector 

v МУ, 

у. Ми, 

W, CW, 

imaginary part of [κ - ω2ν + j ως] 

real part of [к Е ЖЫ! + ісі 

Wy FMW, 

excitation frequency 

time coordinate for graphs 

initial condition of displacement modified to account for 
initial displacement of steady state solution 

initial TDI EE vector 

coefficient matrix of X(0) cos/gfrt term, output of option 1 


coefficient matrix of X(O) sin @,t term, output of option 2 


ordinate of graphical output 


notation lists the principal array and parameter names used in 


the main body of the program. Array names not listed are used for inter- 


mediate operations. Where appropriate the names are associated with the 
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symbols used in the mathematical analysis. 





APPENDIX B 
INSTRUCTIONS FOR USE OF PROGRAM 

B.1 Purpose. The purpose of Program VIBSYS is to determine the natural 
frequencies, mode shapes, and time behavior of a subcritically damped, 
linear dynamic system with N degrees of freedom. (N Z 10) The mass, 
damping and stiffness matrices of the system must be known. 

B.2 Input Data. A blank card must follow the last end card of the pro- 
gram deck. The data cards then follow in the order and format listed be- 
low. The first data card uses the input format 9I5 and ís the control 


card for the program. The nine fields are designated as follows: 


1. М = the order of the system 
22 MPl = 1, Option 1 executed 
0, Omit 
3. | MP2 = 1, Option 2 executed 
0, Omit 
4. MP3 = 1, Option 3 executed 
0, Omit 
5. МР4 = 1, Option 4 executed 
0, Omit 


6. MP5 = 1, Option 5 executed 


0, Omit 

7. ІС - 1, Initial conditions given 
0, Initial conditions are zero 

8. ТЕС = 1, Amplitude and frequency of excitation given 
0, Omit 

9. NPTS Number of points to be calculated for graphical 
output 


All nine fields must be right justified. 





The control card is followed by the mass, damping, and stiffness 
matrices respectively. Each matrix is read in rowwise using input format 
8F10.3* with each row starting on a new card for systems with N Z 8. For 
systems with N > 8 the elements are read in rowwise in a sequential manner, 
that is the first element of the second row should appear in the field ad- 
jacent to the last element of the first row. In this case each matrix must 
start on a new card. 

The initial displacement vector, initial velocity vector, real part 
of the driving force amplitude vector, and imaginary part of the driving 
force amplitude vector, follow, in that order, using format 8F10.3.* 

Each vector must start on a new card, 

The final data card uses a 2F10.3* format. The first field contains 
the excitation frequency and the second the step size (i.e., time increment) 
for the graphical output option. 

Blank cards must be used for initial conditions of displacement and 
velocity, real and imaginary parts of forcing function amplitudes and 
final data card when these values are zero. Multiple runs may be process- 
ed by placing additional data decks behind the first. The program is 


terminated by placing a blank card behind the last data card. 





de 
Should this format present undesirable restrictions on the size of 
the elements of the matrix the so called "E" field may be used by changing 
card number 32. The "E" field was avoided since it is more susceptible to 
error by the user. i 
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B.3 Deck Assembly. The first card of the program deck is a job 
card. The statement "Use scratch tape" must be included in addition to 


the tequired job card information. 


JOB CARD 


PROGRAM VIBSYS 
END 

SUBROUTINE INVERT 
END 


SUBROUTINE MATSUB 
END 
END 


BLANK CARD 


DATA CARDS 


BLANK CARD 
B.4 Cautions to User. The curves drawn in the graphical output option 
are straight line approximations between computed values. The step size 
must therefore be chosen appropriately in order to obtain a smooth curve. 


Since the maximum number of points permitted by subroutine DRAW is 900,it 
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may not be possible to see the transient phase completely die out. 

If difficulty is encountered in obtaining the eigenvalues and eigen- 
vectors, the control parameters of MATSUB may be altered to enhance con- 
vergence. The value of MIT and MITS, card numbers 20 and 21, control the 
maximum number of iterations to be performed in the power method and in- 
verse power method respectively. ALRS and ALIS, card numbers 78 and 79, 
represent a complex parameter that may be called the "origin". MATSUB 
will usually converge on the eigenvalue most distant from the "origin". 
ALRS = 1.0 and ALIS = 0.1 in Program VIBSYS. The convergence criteria 
for the power method is set by the value of EPl, card number 85 while 
the inverse power method is determined by EP2 card number 86. ЕРІ - 107°, 
and EP2 = 10714 in the present version. Furthermore, if IEG, card number 
76, is set equal to one, the eigenvalue iterants will be printed out and 


more suitable values of ALRS and ALIS may possibly be found by inspection. 
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APPENDIX C 


FLOW CHARTS 


The symbols used in the flow charts are defined as follows: 


- arithmetic, read or print operations 


: connectors 


- decísion 


- call subroutine indicated 


0991 


The notation used in the flow chart for MATSUB corresponds to that used 


in the CO OP Library version, reference 7. 





FLOW CHART MAIN EODY PROGRAM VIBSYS 


READ: PROGRAM PARAMETERS, N, MP1, MP2, 





READ: MASS, DAMPING, AND 
. STIFFNESS MATRICES 


ASS, DAMPING, AND 


πι 
STIFFNESS MATRICES 


TREAD: INITIAL CONDITIONS | 
XO, VO 





PRINT: XO, VO 


READ: DRIVING FORCE AMPL. 
DR, DI 


PRE OR, DI 


READ: EXCIT. FREQ. AND 
STEP SIZE. WO, STEP 


PRINT: WO, STEP 


READ: ITERATION PARAMETERS 


ND OF INPUT DATA 








AMM- AM . 
LES | 
E 14 
AMS- —AM*S 


AMD = - AM*DM 


FORM UR MATRIX 


CALL MATSUB 
UR, UI 


PRINT: TRACE, DET, EIGENVALUES, 





EIGENVECTORS, SUM, AND PROD 


\ 
SELECT POSITIVE VALUES OF VALI 
[SPEC = VALI**2- VALR**2 ——— | 


31 


WNAT = SQRTF( |I| SPECII ) 


SELECT VECR, VECI CORRESPONDING TO 
POSITIVE VALI 


DISPR = VECR 
DISPI = VECI 


PRINT: NATURAL FREQUENCIES 
AND MODE SHAPES 








940 
YES NO 
951 
PRINT WNAT = O SYSTEM IS DIF = WNAT(I)- WNAT(J) 
: OVERDAMPED | 
= NO YES 
START FORCED VIBRATION WITH | 941 
SINUSOIDAL EXCITATION PRIN Y Y HAS 


O AL Mach 


WIR = S - WO**2*AMM GO TO 1000 
wer 0 

WDM =WO* DM 

WMKI = ММК 


ALAR 


WIC = T m 


904 


807 
= Y A | s 3 


А [4] 


910 
RMI = RTD ~ ТІРІ 
CPD = RTD + TTD 


L i2 UN PROBLEM 





I=1,N 
VMV = DISPR*AMM*DISPR 
WMW = DISPI*AMM*DISPI 
VMW = DISPR*AMM*DISPI 
WCW = DISPI*DM*DISPI 


VCV = DISPR*DM*DISPR 
VCW = DISPR*DM*DISPI 
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GR = B/ARN**2 * BRN**2 
HR = -HM/ARN**2 * BRN**2 





81 
CRM - 2* G3 * AMM 
HRM = 2*HR*AMM 
GRC = 2*CR*DM 
HRC = 2*HR*DM 
82 


ХОС = = AGM + EHM + GRC 
XOS = -BGM - AHM + HRC 


NO YES 
101 









Me 


201 
АВ = GRM*VO + XOC*XO 
CD » HRM*VO * XOS*XO PRINT  XOC,XOS 
. | co TO 100 


YES MP] = 0 
KC = 2 . Ä 
CO TO 60 


PRINT  AB,CD 
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ACO - XO - RMI 
VCO = VO = WO*CPD 









710 
- GRM*VCO * XCC*XCO 
CDC = HRM*VCO + XOS*XCO 





711 


NO | σετ YES 














CALCULATE 
RANS TER 





GALCULATE GRAPH POINTS FOR 


SIEADY STATE SOLUTION <> 


FIND MAX ORDINATE 
TEST MAX ORDINATE 


Я 719 
1 Te 
"n 







a 


YSCALE 







MAX ORD/2 














STORE GRAPH TITLES 
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SUBROUTINE MATSUB 


SET: M, IEG, IVEC, «, , 
О, DUE MISC C 


UH = C 
PROD = 0 
COMPUTE TRACE O 


ASSIGN 520 TO IA 
ASSIGN 811TO 1D 
INTER 0 


TRIANGULARIZE A 


COMPUTE DET À 


PRINT TRACER,TRACEI , 
DET A 








ASSIGN 912 TO ID 
ASSIGN 530 TO IA | 
ASSIGN 40 TO IB 


AN 





© BA 
JSL < 0 
Q 
[17:1] 
I3 1.0 
Q А-е 1) > А 
1-1 





Е 


^ 


ITERATION POWER METHOD 
L = (A - e 1 Pe 


= E oro YES 
gk = 1000 


ш - NO 
JAS от ш | 
471) „el 
| Ἡ ERROR PRINT 





— EET 


PRINT 1, , 44; 


PRINT N,ALR,ALI,ICT,IJ 
y 


VALR(N) = ALR 
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YES 
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SET FLAG 
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NO 
П ХЦ 


214 
ISL = 0 


01 
NO YES 






| VALR(N) = ALR 
VALI (N) = ALI 






SAVE &i 
SUM = SUM + ©: LOC 
DROD*# oti 65 | PRINT NOTE AND IT č | 
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с Y 990 

. | MATRIX DEFLATION É ["" 40: EXIT 
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B=DEFLATED A 
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SUM = SUM + <; 


PRINT SUM, PROD Í PROD* % 
ος 
5 Б 990: 











935 


TRIANGULARIZE A 
| 911 


TIEN = 1 | 
ASSIGN 530 TO IA 
920 


COMPUTE DET A 
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PRINT DET A 


924 


COMPUTES AND PRINT LARGEST 
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4 е 
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BACKWARD SUBSTITUTIO 
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PRINT EIGENVECTOR 
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APPENDIX D 


PROGRAM LISTING 
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APPENDIX E 
SAMPLE PROBLEM 
The sample problem presented is eut of a four degree of freedom 
system shown in the figure below. In order to show a comparison of data 
for damped and undamped systems, the example chosen is one in which the 
natural frequencies of the undamped system may be obtained with relative 


ease. The program is not restricted to problems with the symmetry dis- 


played in the example. 


f) 
| X, ye Ka hos mo^ 


ib F: Bb E 


2 
10 1b sec/in 


Let HM I 
C1 a C^ T C4 ы σι = 30 lb sec/in 


For the undamped case the frequency determinate then becomes 


2K - Mw? -K 0 0 
-K ж - мо? -K 0 
0 -K 2K - MW? -K 
0 0 -K 2K - Mw* 





Expanding the determinate 


8 K Ё K : 4- < 5 2. K il 

4% - в) + ait) eo - 200) бу + s = © 
The roots of the frequency equation are: 
Q, = 6.18034, 02, - 11.75571,45: - 16.18034, 024 - 19.02113 

The problem was programmed with a forcing function f(t) = 100віп 
(60t), and an initial displacement on M,of 0.5 inches. A step size of 
0.005 seconds was used for the graphical option. All five options were 
called for and the resulting output is shown on the following pages. 

In addition to the output presented for C = 30 lbsec/in runs were made 
with C = 0.01, 0.10, and 1.0 1b sec/in. The frequencies and damping 


ratios for all runs are tabulated below. 








C 1b sec/in 5 (9, , rad/sec ЕСА C2, ,rad/sec 
0.01 — 6.18034 0.0001 11.75571 
0.10 0.0003 6.18034 0.0006 11.75570 
1.00 0.0031 6.18028 0.0067 11.75530 
30.00 0.0098 6.12699 0.1821 11.38430 

C 1b sec/in Be ر وا‎ rad/sec ω, rad/sec 
0.01 0.0001 16.18034 0.0001 19.02113 
0.10 0.0008 16.18033 0.0009 19.02111 
1.00 0.0081 16.17928 0.0095 19.01941 
30.00 0.2584 15.19737 0.3118 17.40396 
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